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Abstract 

We report on a study of QCD thermodynamics with three flavors of quarks, using a Symanzik 
improved gauge action and the Asqtad O(a^) improved staggered quark action. Simulations were 
carried out with lattice spacings l/4r, 1/6T and 1/8T both for three degenerate quarks with 
masses less than or equal to the strange quark mass, iris, and for degenerate up and down quarks 
with masses in the range 0.1ms < m^^d < 0.6ms, and the strange quark mass flxed near its physical 
value. We present results for standard thermodynamics quantities, such as the Polyakov loop, the 
chiral order parameter and its susceptibility. For the quark masses studied to date we find a rapid 
crossover rather than a bona fide phase transition. We have carried out the first calculations of 
quark number susceptibilities with three fiavors of sea quarks. These quantities are of physical 
interest because they are related to event-by-event fluctuations in heavy ion collision experiments. 
Comparison of susceptibilities at different lattice spacings show that our results are close to the 
continuum values. 

PACS numbers: 12.38.Gc, ll.lO.Wx, 12.38.Mh 
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I. INTRODUCTION 



Under ordinary laboratory conditions quarks and gluons, the fundamental constituents 
of quantum chromo dynamics (QCD), are bound into hadrons. However, as the temperature 
or density is increased, the forces among quarks and gluons weaken, and one expects to find 
a phase transition or crossover to a new state of matter, a quark-gluon plasma. The plasma 
was a dominant state of matter in the early development of the universe, and a primary 
objective of current relativistic heavy ion collision experiments is to observe the plasma and 
determine its properties. The behavior of strongly interacting matter in the vicinity of the 
phase transition or crossover is inherently a strong coupling problem, which at present can 
only be addressed from first principles through lattice gauge theory calculations. Among the 
issues that can uniquely be addressed by lattice calculations are the nature of the transition 
and the temperature at which it occurs, the properties of the plasma, including strange 
quark content, and the equation of state 

In this paper we report on a study of the phase diagram of high temperature QCD with 
three flavors of quarks using improved gauge and quark actions Most lattice studies of 
high temperature QCD have included only the up and down quarks, but the inclusion of 

n 

the strange quark is an important feature of our work 3]. From chiral models it is expected 
that with two flavors of quarks there is no phase transition at all for physical values of the 
up and down quark masses; however, a strange quark could induce a flrst order transition, 
or move a second order critical point closer to physical quark masses ^. Furthermore, the 
production of strangeness is expected to be an important signal for the plasma in heavy ion 
experiments. 

We are considering two cases in our study: 1) all three quarks have the same mass m^; 
and 2) the two lightest mass quarks have equal mass m^ ^, while the mass of the third 
quark is flxed at approximately that of the strange quark m^. We refer to these cases as 
Nf = 3 and Nf = 2 + 1, respectively. For Nf = 3 we have worked with quark masses in 
the range 0.2 ms < ruq < rus, while for Nf = 2 + 1 we have carried out simulations with 
0.1 < rriu^d < 0.6 m^. We have monitored standard thermodynamic quantities such as 
the plaquette, the Polyakov loop, the chiral order parameter and its susceptibility. For the 
quark masses we have studied to date, we flnd rapid crossovers, which sharpen as the quark 
mass is reduced, rather than a bona fide phase transition. This is in agreement with earlier 
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work We have also measured the quark number susceptibihties, which provide excellent 
signals for the crossover, and which are directly related to event-by-event fluctuations in 
heavy ion collisions 0]. The data indicate that our results for these susceptibilities are quite 
close to their continuum values. 



II. THE SIMULATIONS 

Our simulations are carried out with a one-loop Symanzik improved gauge action |^ and 
the Asqtad quark action [3]. Both the gauge and quark actions have all lattice artifacts 
removed through order a} at tree level, where a is the lattice spacing, and are tadpole 
improved. Thus, the leading order flnite lattice spacing artifacts are of order a^a^ and a^. 
We and others have used this action over the last few years to perform a wide variety of 
zero-temperature calculations ^. It has a number of features that make it particularly well 
suited for high temperature studies. The Asqtad action has considerably better dispersion 
relations for free quarks than the standard Kogut-Susskind and Wilson actions which 
markedly decreases lattice artifacts above the transition. This is also jtrue of the P4 action 
used by the Bielefeld group in its studies of high temperature QCD IQj. The improvement 
in dispersion relations is illustrated in Fig. ^ where we plot the energy, pressure and quark 
number susceptibility for free massless quarks as a function of lattice spacing. The Asqtad 
action exhibits excellent scaling properties in the lattice spacing 11], which accelerates the 
approach to the continuum limit. Finally, taste symmetry breaking is much smaller for 
the Asqtad action, than for the conventional Kogut-Susskind action. Our spectrum studies 
with the Asqtad action indicate that for lattices with eight to ten time slices, the kaon is 
heavier than the heaviest non-Goldstone pion in the neighborhood of the flnite temperature 
transition or crossover, a condition which requires much smaller lattice spacings with the 
conventional Kogut-Susskind action. It is, of course, difficult to study the effects of the 
strange quark on the transition if this condition is not fulfilled. 

We have used the refreshed hybrid molecular dynamics R algorithm to generate gauge 
configurations. In integrating the molecular dynamics equations of motion we used a time 
step dt equal to the smaller of 0.02 and \mu^d- The momenta conjugate to the gauge fields 
were refreshed every molecular dynamics time unit, which consisted of Xjdt time steps. In 
the vicinity of the crossover we generated 2,000 equilibrated molecular dynamics time units 
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at each value of the gauge couphng and quark masses we studied. The runs well above and 
below the crossover were shorter. Measurements of standard thermodynamics quantities 
were made after each time unit, and gauge configurations were saved every five time units 
for separate measurements of the quark number and iptp susceptibilities. Ten random sources 
were used in the measurements of ipip niade at the end of each time unit, and 100 random 
sources were used in the calculations of the susceptibilities. The values of the chiral order 
parameter obtained in these two measurements agreed to within statistics. 

We have attempted to vary the temperature while keeping all other physical quantities 
constant. To this end, for the Nf = 3 study we have performed a set of spectrum calculations 
at lattice spacings a = 0.125 fm and 0.18 fm with nig = rris, 0.6 m^, 0.4ms, and 0.2 ms. 
We determine the lattice spacing from the static QQ potential, and express it and other 
dimensionful quantities in terms of ri defined by tIFqqIti) = 1. Using results from the 
1S-2S and IP-IS splitting in the T spectrum [l^, we have found ri = 0.317±0.007 fm in the 
continuum and chiral limits Q|. We determine m^ from the requirement that m^^Jm^ ~ 
0.673, where rjss is an "unmixed" pseudoscalar meson made of an s and s quark. We would 
like to carry out our thermodynamics studies with three equal mass quarks for m^^^/m^ 
fixed, but this quantity will, of course, vary slightly with lattice spacing if we keep mg/m^ 
fixed. So, at a = 0.18 fm we perform linear interpolations of m^^^ and m^ in the quark mass 
to determine the precise values of for which mr^^Jm^ will take on the values found at 
a = 0.125 fm. Then to determine the values of arUg, a and T = l/aNf for thermodynamics 
studies with other values of a, we perform interpolations or extrapolations with a form due 
to AUton, inspired by asymptotic freedom jl^. We first define 



where 6o aiid bi are the universal beta-function coefficients for massless three-flavor QCD, 
and is the bare lattice coupling, fig"^) is basically q/Al. We then determine a as a 
function of g^ from the interpolation formula 



The second term in Eq. Q is an 0{g^ o?) correction to the asymptotic freedom formula. 
The coefficients Cq and C2 are determined from the measured values of a/r i at a = 0.18 and 



(1) 




(2) 
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0.125 fm. A similar interpolation formula is used to determine the lattice quark mass along 
the line of fixed m^^^/m^, 

am,{g') = do {bogY'^' fig') [l + d,g'f{g')] . (3) 

Here we have included the anomalous dimension of the mass. In the initial stages of our 
work we simply made linear interpolations of ln(a) and In(amq) in l/g"^ between the anchor 
points at a = 0.18 and 0.125 fm, which gave results in good agreement with those of Eqs. (j21) 
and 0. However, use of the more sophisticated formulae became critical for points outside 
the anchors. 

Our approach for thermodynamics studies with rriu^d < n^'s is quite similar. In this case 
we wish to vary the temperature keeping both niT^/mp and m^^Jm^ fixed. We have carried 
out spectrum studies at a = 0.125 and 0.18 fm for light quark masses mu,d = f^s, 0.6 m^, 
0.4 ms, 0.2 ms, and 0.1 m^, and at a = 0.09 fm with m,, d = m^, 0.4 ms 0.2 ms and 0.1ms 
all with the mass of the heavy quark fixed close to rris jsl, |lj| . The strange quark mass was 
determined from the spectrum calculations with three equal mass quarks. In our spectrum 
runs at a = 0.125 fm 9| we found that m^^^ varied by less than 2% and m,^ by less than 
1% for 0.2 ms < m„_rf < rus, and the heavy quark mass held fixed at rris. So, the neglect 
of the dependence of rrir,^^ and m^ on mu,d is well justified. In these studies we performed 
linear interpolations of m^ and rup at a = 0.18 fm to determine the values of m^ ^ for which 
mTj/nip takes on the values found at a = 0.125 fm. Then, for other values of a we interpolate 
or extrapolate a/ri using Eq. and arriu^ and anis using Eq. (jH)). We again determine 
the values of Cq, C2, do and d2 from the spectrum runs at a = 0.125 fm and 0.18 fm. For 
TTtu,d = ^s, 0.4ms, 0.2 ms and 0.1 rris, where we have spectrum data at a = 0.09 fm, we can 
add f^ig') terms to the right hand sides of Eqs. ^ and ©. These added terms make only a 
few percent difference in the lattice spacing and quark masses extrapolated to lOT,^, giving 
us confidence in the interpolations and extrapolations used in our study. We have recently 
combined data from all zero temperature runs we have made to date to obtain a smooth 
interpolation formula of ln(ri/a) as a function of quark mass and gauge coupling The 
results are in excellent agreement with those obtained from Eq. (2). 
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III. RESULTS FOR Nf = 3 



For three equal mass quarks, Nf = 3, we have carried out thermodynamics studies on 
lattices with four, six and eight time slices, and aspect ratio Ng/Nf = 2. Here A^^ and 

are the spatial and temporal dimensions of the lattice in units of the lattice spacing. 
We also performed simulations with aspect ratio three for A^^ = 4, and obtained results 
that are indistinguishable from those with aspect ratio two. The spectrum calculations and 
interpolations described above allowed us to determine the values of the quark mass, ruq, 
that keep m^^^ / fixed as the gauge coupling is varied. They also enabled us to determine 
the value of the lattice spacing, and therefore the temperature for each run. 

In Fig. 121 we plot the real part of the Polyakov loop as a function of temperature on 16'^ x 8 
lattices. In this figure, as elsewhere, we give the values of mq/rris for a = 0.125 fm. The 
corresponding values of m^^^/m^ are given in Table lU The Polyakov loop shows a crossover 
from confined behavior at low temperature to deconfined behavior at high temperature. 
There is a slight trend for the temperature dependence of the Polyakov loop to be steeper 
for larger quark masses. This is to be expected, since at sufficiently large quark masses, it 
is a bona fide order parameter. The insensitivity of the Polyakov loop to the quark mass 
at fixed temperature or lattice spacing is, perhaps, not surprising. We have determined the 
lattice spacing from the heavy quark potential, and in our spectrum runs we adjusted the 
coupling constant to keep the lattice spacing fixed as the quark mass is varied. Since the 
Polyakov loop, like the heavy quark potential, is determined from measurements of purely 
gluonic operators, our procedure is likely to minimize the dependence of the Polyakov loop 
on the quark mass. 



The number of conjugate gradient iterations varies as the temperature is changed refiect- 
ing a changing condition number for the Dirac operator as the spectrum of states changes. 
We show this quantity in Fig. IHl for the four values of the quark mass we have studied on 
16^ X 8 lattices. The sharpening of the crossover as the quark mass is reduced is evident. 
The ipip susceptibility provides a more physical signal for the crossover. It is given by 

Xtot = = tt{Tt M-'Tt M"1) - -(Tr M'^) - -(Tr M^Y , (4) 

am Vs Vs Vs 

where M is the fermion matrix, T the temperature, and Vg the spatial volume. The traces 
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= 3 


Nj = 


2 + 1 


rriq/ms 




mu,d/ms 




1.0 


0.673 


1.0 


0.673 


0.6 


0.583 


0.6 


0.582 


0.4 


0.504 


0.4 


0.509 


0.2 


0.404 


0.2 


0.392 






0.1 


0.298 



TABLE I: In the first column we show the vakie of mq/ms at lattice spacing 0.125 fm, which 
produced the "rnj^^^/m^ ratio shown in the second column for spectrum calculations with three 
equal mass quarks. In the third column we give the value of rriu^d/'nis which produced the m-^jmp 
ratio shown in the fourth column for spectrum calculations with two equal mass light quarks, and 
the mass of the heavy quark fixed at m<j. 

in Eq. (jH) were, as usual, evaluated using the identity 

Tr O = (i?*Oi?)H, (5) 

where Or indicates an average over vectors, i?, of Gaussian random numbers. We used 100 
vectors of random numbers for each gauge configuration. Of course, care must be taken not 
to use the same vector of random numbers in evaluating the product of traces in the first 
term of Eq. (0}. We plot Xtot in Fig- El for quark masses m^, 0.6 ms, 0.4 ms, and 0.2 ms on 
12^ X 6 lattices. Note the increase in the height of the peak as the quark mass is decreased. 

In Figs. El and El we show the chiral order parameter, ?/'?/' as a function of temperature 
for vTiq = rris, 0.6 ms, 0.4 ms, and 0.2 ms on 12^ x 6 and 16'^ x 8 lattices, respectively. The 
bursts in these figures are linear extrapolations of ipip in the quark mass to rriq = for fixed 
temperature. These figures suggest that for rriq = there is unlikely to be a phase transition 
for temperatures above 190 MeV, but one could occur at or below that value. For sufficiently 
high temperatures and small quark masses, one expects a linear extrapolation in the quark 
mass to be valid, so the vanishing of ipip at rriq = indicates that the system is in the chiral 
symmetric phase. However, a non-vanishing value of the extrapolated tpip could indicate a 
first order transition at a non-zero quark mass, or simply a breakdown in the validity of the 
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linear extrapolation in the quark mass. 



IV. RESULTS FOR Nj = 2 + 1 

The Nf = 2 + 1 thermodynamics studies were carried out primarily on 8^ x 4, 12^ x 6 and 
16^ X 8 lattices. In addition, a number of runs were made on 18^ x 6 lattices to check for finite 
spatial size effects, which turned out to be negligible. In this phase of our work, we performed 
simulations with two degenerate light quarks and the heavy quark mass held approximately 
equal to that of the strange quark. The spectrum calculations and interpolations described 
above allowed us to determine curves of constant physics in the coupling constant-quark 
mass plane, and to determine the lattice spacing, and therefore the temperature for each 
run. 

In Fig.[7|we plot the real part of the Polyakov loop as a function of temperature on 16'^ x 8 
lattices for the four values of m^^rf studied to date. As in the Nf = 3 study, we observe a 
crossover from confined to deconfined behavior, rather than a sharp transition, and little 
dependence on the light quark mass. In Fig. |Hlwe show the number of conjugate gradient 
iterations needed for inversion of the Dirac operator for the light quarks on 16'^ x 8 lattices, 
and in Fig. 1^1 we plot the ipip susceptibility for the up and down quarks on 12'^ x 6 lattices. In 
both figures the sharpening of the crossover as the light quark masses are reduced is evident. 

In Figs. El and ^2 we plot the chiral order parameter as a function of temperature on 
12^ X 6 and 16^ x 8 lattices. The octagons, diamonds, squares and fancy crosses are data for 
= 0.6 rris, 0.4 m^, 0.2 and 0.1 respectively, and the bursts are linear extrapolations 
in rriu^d for fixed temperatures using the 0.2 ms and 0.1 rris points. The fact that the linear 
extrapolation of the chiral order parameter becomes slightly negative suggests that in this 
region of temperature and quark mass ipip is nonlinear, as would be expected if there were 
a bona fide critical point for rriu^d > 0. 



V. QUARK NUMBER SUSCEPTIBILITIES 



In order to study the quark number susceptibihties 



16 



17 



19 



20 



2l|, we introduce 



chemical potentials coupled to a set of mutually commuting conserved charges Qa- The 
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partition function can then be written 



exp 



(6) 



The quark number susceptibihties can be related to event-by-event fluctuations in heavy ion 
collisions Q by the fluctuation-dissipation theorem 



(7) 



Vs d Had fx (3 

We work at /i^ = 0, so the brackets, (), in Eq. ((7j) indicate averages weighted by the standard, 
real Euclidean action for QCD, and (Qa) = 0. The charges we consider are bilinears in the 
quark fields, so it is possible to write Xo,/3 in the form 



Xa,l3 — Ca^p + D, 



a,l3 



(8) 



where C^,/? is the contribution to Xo,/3 of connected diagrams in which quark fields from 
charge Qa contract with those from Q^, whereas -D^,/? is the contribution of the disconnected 
diagrams in which quark fields within each charge contract. 

If we take the charges to be the number operators for up, down and strange quarks, then 
because we have set m„ = m^, x takes the form 



X 



Ci^i + Di^i Di^i 
Dii Cii + Dii 



\ 



\ 



Dig Di s Cs s + Dg s 



(9) 



where the rows and columns of the matrix are labeled by m, d and s in that order, and / 
stands for the equal up and down quark matrix elements. 

It appears more physical to take the three independent charges to be the z-component of 
isospin, Q/, the hypercharge, Qy, and the baryon number, Qb, where in the continuum 



Qi 

Qy 
Qb 



x] 



V3 



(10) 
(11) 
(12) 



10 



Here is a three-component column vector, whose components are the up, down and 
strange quark fields, and A3 and As are the standard diagonal generators of SU(3) in the 
fundamental representation, A3 = diag(l,— 1,0) and Ag = diag(l, 1, — 2)/-\/3. Using these 
charges, x can be written in the form 



^ \Ci,i ^ 

x= |(Q,z + 2C,,,) + |(A,^-2A,s + /^s,s) i(a,/-C,,, + 2A,z- A,s-/^.,s) 

^ l{Ci,i-Cs,s + 2D^-Di^s-Ds,s) \{2Ci,i + Cs,s + m,i+m,s + Ds,s) ) 

(13) 

where the rows and columns are now labeled by /, Y and B in that order [22]. Note 
that in the Nf = 3 simulations for which = = m^, Ci^i = Ci^s = Cs,s and Di^ = 
Di s = Dg s, so X is a diagonal matrix in this representation. There are then no correlations 
between fiuctuations /, Y and B. In the Nf = 2 + 1 case the only correlations are between 
hypercharge and baryon number. Of course for m„ 7^ there would be correlations 
among the fiuctuations in all three charges. For temperatures below the phase transition 
or crossover, the lightest particle that can be excited by a chemical potential coupled to 
the z-component of isospin is the pion, while for hypercharge and baryon number chemical 
potentials it is the kaon and the nucleon, respectively. Above the transition temperature 
each of the chemical potentials can excite quark states that are much lighter than hadrons, 
so we expect the diagonal elements of x to increase sharply in the vicinity of the transition, 
and they do. 

As in the case of Xtot, these susceptibilities can be written as expectation values of traces 
of the quark matrices and their derivatives with respect to the chemical potential, and 
unbiased estimators of the traces can be expressed in terms of vectors of Gaussian random 
numbers Q]. Here too, we use 100 random vectors for each gauge configuration. 

In Fig.^we plot the quark number susceptibility of massless free quarks for the standard 
Kogut-Susskind, Wilson, P4 and Asqtad actions. (With only one fiavor of quark, there is, of 
course, only one susceptibility). One sees that as in the case of the energy and pressure, the 
susceptibility of the Asqtad action is significantly closer to the continuum result for small 
values of Nf (large values of the lattice spacing) than that of the Kogut-Susskind action. 
In Fig. ^1 we show xi,i for the four quark masses we have studied on 12'^ x 6 lattices with 
Nf = 3. One again sees the steepening of the crossover and its shift to lower temperature 
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as the quark mass is decreased. In Fig. [121 we illustrate the dependence of xi,i on lattice 
spacing by plotting results at fixed quark mass for three values of A^*. The solid lines on the 
right of this figure indicate the values of the susceptibility for free quarks in the continuum 
and on the finite lattices on which the simulations were carried out. 

In Fig. [131 we plot xi,i on 12^ x 6 lattices for Nf = 2 + 1. As indicated in Eq. (fT^. 
this quantity is proportional to the light quark connected diagram, Ci^i. The disconnected 
diagrams are, of course noisier. As an example, we plot the light quark disconnected graph 
in Fig. El This quantity is given by 

A,^ = {{B + Y/2Y) - {{B + Y/2)f - xi,i. (14) 

We see in Fig. El that Di^i can be cleanly evaluated in the neighborhood of the crossover, 
the only region in which it is appreciable. Below the crossover one expects xi,i to be larger 
than Xb+y/2,b+y/2i since in this regime the lowest energy state that can be excited by a 
chemical potential coupled to is a pion, while the lowest energy state that can be excited 
by a chemical potential coupled to S + F/2 is a kaon. Our results suggest that vestiges of 
hadronic physics persist in the plasma at least up to 240 MeV. In Figs. El and El we show 
results for xi,i and for the strange quark number susceptibility, Cs,s + -D^^s, on lattices with 
4, 6 and 8 time slices for a range of spatial volumes. There is no observable dependence 
on the spatial volume. The close agreement between the Nf = 6 and 8 results here and 
for Nf = 3 illustrates the excellent scaling properties of the action, and indicates that our 
results are close to the continuum ones. Finally, in Fig. El we plot the diagonal elements 
of the susceptibility matrix, xi,iy Xy,y and Xb,b, as a function of temperature for two light 
quarks with mass 0.2 rris and one heavy quark with mass on 12^ x 6 lattices. Xy.y and 
Xb,b have been multiplied by factors of 3/4 and 3/2 respectively, so that the quantities 
plotted approach the same high temperature limit as xi,i- Also shown is Xy,b, the only 
non-zero off-diagonal matrix element of x for rriu = md- It measures correlations between 
fluctuations in the hypercharge and baryon number. The coefficient of xi,b in this figure is 
the geometric mean of those for xy,y and Xb,b- 
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VI. CONCLUSION 



For the quark masses we have studied to date, we find a rapid crossover, rather than 
a bona fide phase transition, for both Nf = 3 and Nf = 2 + 1. Our resuh for Nj = 3 is 



consistent with recent work of Karsch et al. 



23| . They find a first order phase transition for 



Nf = 3 only for pion masses below 290(20) MeV for the standard Kogut-Susskind action 
and below 67(18) MeV for the improved P4 action on lattices with four time slices. The 
lightest quark mass used in 



23| corresponded to a pion mass of 170 MeV, and reweighting 
techniques were used to extrapolate to lighter quark masses. Our lightest Nf = 3 quark 
mass corresponds to a pion mass of approximately 340 MeV. Because of the large difference 
between standard and improved actions on these lattices, it seems particularly important to 
push our work at Nt = 6 and 8 to smaller quark masses, and we intend to do so. 

The small quark mass that seems to be required for a first order transition at A'^^ = 3 
strongly suggests that in the real world, Nf = 2 + 1, there is no phase transition at the 
physical quark masses. With this assumption, we have estimated the critical temperature 
for Nf = 2 + 1 at m„ = through an extrapolation of the form 



ri Tc = Co + ci {m^/mpY + C2 (aTj^, (15) 

where we evaluated for each value of Nt and m^/ m^, for which we have made measurements 
from the peak in the iptp susceptibility. For a second order phase transition in the 0(4) 
universality class at mu,d = 0, d = 2/135 ^ 1.08. We find that = 169(12)(4) MeV with 
a of 2.1 for 11 degrees of freedom. The first error is the fit error, the second from the 
uncertainty in ri, taken as 0.317(7) fm jlj]. To test the sensitivity of Tc to d, we have also 
performed a fit with d = 2, which yields Tc = 174(11) (4) MeV with a of 1.5 for 11 degrees 
of freedom. So, the goodness of the fit does not allow us to prefer either of them. 

Finally, we note that the agreement between the quark number susceptibilities on Nf = 6 
and 8 lattices is very encouraging, as it indicates that our results for these quantities are 
close to their continuum values. 
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FIG. 5: The chiral order parameter, (ipip), as a function of temperature on 12^ x 6 lattices for 
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FIG. 7: The real part of the Polyakov loop as a function of temperature for two light and one 
heavy quark on 16^ x 8 lattices. Results are shown for light quark masses rriu^d/fns = 0.6, 0.4, 0.2 
and 0.1. The mass of the heavy quark is fixed at m<j. 
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FIG. 8: The number of conjugate gradient iterations required for convergence of the inversion of 
the Dirac operator of light quarks on 16'^ x 8 lattices. Results are shown for light quark masses 
"i«,d/"is =, 0.4, 0.2 and 0.1. The mass of the heavy quark is fixed at m^. 
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FIG. 9: The ipip susceptibility as a function of temperature for two light and one heavy quark on 
12^ X 6 lattices. Results are shown for light quark masses mu^d/rns = 0.6, 0.4, 0.2 and 0.1. The 
mass of the heavy quark is fixed at m^. 
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FIG. 10: The chiral order parameter ^^/^ as a function of temperature on 12^ X 6 lattices. The 
bursts are linear extrapolations in the quark mass to rriud = for fixed temperature. 



25 



A 
V 

CO 



0.15 



N.=2+l, le^'xS 



0.10 



0.05 - 



0.00 







□ 



o 



© 

□ o ® 

X m ^ o 

□ ^ o 



"1 — I — r 



° m^,d=0. 61113 
o m^d=0-4m3 



O 



SIX 

✓ I \ 



X 



\ I ✓ 
-I- 







X 



-I- 



□ 







\ I ✓ 
✓ I \ 



o 



^1^ /l\ 



J L 



XX X X X X X y. 



J I L 



J I I L 



J I L 



150 



200 250 
T (MeV) 



300 



FIG. 11: The chiral order parameter tpxp on 16^ x 8 lattices. The bursts are linear extrapolations 
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FIG. 12: The z-component of isospin susceptibility xi,l ^ ^ function of temperature for three 
degenerate flavors of quarks with masses mq/ms = 1.0, 0.6, 0.4, and 0.2 on 12'^ x 6 lattices. The 
solid lines on the right of the figure indicate the free quark value in the continuum and on a 12^ X 6 
lattice. 
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FIG. 13: The z-component of isospin susceptibility xi,i for three degenerate flavors of quarks with 
mass nig = 0.4ms, on 8'^ x 4, 12^ x 6 and 16'^ x 8 lattices. The sohd hues on the right of the 
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simulations were carried out. The close agreement among these results is another indication of the 
excellent scaling properties of the Asqtad action. 
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FIG. 14: The z-component of isospin susceptibility as a function of temperature for two light and 
one heavy quark on 12'^ x 6 lattices. Results are shown for light quark masses mu^/ms = 0.6, 0.4, 
0.2 and 0.1. The mass of the heavy quark is fixed at m^. The solid lines on the right of the figure 
indicate the free quark value in the continuum and on a 12^ x 6 lattice. 
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FIG. 18: The diagonal elements of the susceptibility matrix, x/,/) XB,B and XY,Y as a function of 
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